function [modelYield] = projectionStepBondPricesClosedForm(model,setup)

maxMat = setup.maxMat;

% Save output for the core of the model to compute bond prices
thetaTrans = model.theta';
setup.out  = saveStatesControls(thetaTrans(:),setup);

% Computing bond prices 
nx = setup.nx;
if setup.appOrder == 1  
    thetaBondPrice      = zeros(maxMat,1+nx);
elseif setup.appOrder == 2
    thetaBondPrice      = zeros(maxMat,1+nx+nx*(nx+1)/2);
elseif setup.appOrder == 3
    thetaBondPrice      = zeros(maxMat,1+nx+nx*(nx+1)/2+nx*(nx+1)*(nx+2)/6);
end
MAEbonds   = zeros(maxMat,1);
RMSEbonds  = zeros(maxMat,1);
MaxAEbonds = zeros(maxMat,1);
for k=1:maxMat
    if k == 1
        setup.out.p1_cup = zeros(size(setup.xGrid,1),setup.n_nodes); %log-transformation so exp(0) = 1
    else
        setup.out.p1_cup = p2_cup;
    end
    % Optimization
    setup.computeP2_cup = 0;
    [~,~,~,logP2_cu] = NLSobjectFunc(thetaBondPrice(k,:)',setup);
    thetaBondPrice(k,:) = OLSforProj(logP2_cu,setup.xGrid,setup.appOrder,setup.scaleX_in_g);
    
    % Computing p2_cup and the residuals
    setup.computeP2_cup = 1;
    [~,resEuler,p2_cup] = NLSobjectFunc(thetaBondPrice(k,:)',setup);
    MAEbonds(k,1)  = mean(abs(resEuler(:)));
    RMSEbonds(k,1) = sqrt(mean(resEuler(:).^2));
    MaxAEbonds(k,1)= max(abs(resEuler(:)));
end

% Computing bond yields (in "by" - for bond yields)
scaleX_in_g    = setup.scaleX_in_g;
nn             = -1./(1:maxMat)';
bySS           = repmat(log(model.params.Rss),maxMat,1);
by0            = nn.*thetaBondPrice(:,1);
byx            = nn.*thetaBondPrice(:,1+1:1+nx)./repmat(scaleX_in_g',maxMat,1);
byxx           = zeros(maxMat,nx,nx);
byxxx          = zeros(maxMat,nx,nx,nx);
if setup.appOrder > 1
    idx   = 0;
    for i=1:nx
        for j=i:nx
            idx = idx + 1;
            byxx(:,i,j)    = nn.*thetaBondPrice(:,1+nx+idx)/(scaleX_in_g(i,1)*scaleX_in_g(j,1));
            byxx(:,j,i)    = byxx(:,i,j);
        end
    end
end
if setup.appOrder > 2
    idx   = 0;
    for alfa1=1:nx
        for alfa2=alfa1:nx
            for alfa3=alfa2:nx
                idx = idx + 1;
                byxxx(:,alfa1,alfa2,alfa3)    = nn.*thetaBondPrice(:,1+nx+nx*(nx+1)/2+idx)/(scaleX_in_g(alfa1,1)*scaleX_in_g(alfa2,1)*scaleX_in_g(alfa3,1));
                
                % Using symmetry for alfa1 and alfa2
                if alfa1 == alfa2 && alfa2 ~= alfa3 %alfa1==alfa2~=alfa3
                    %byxxx(:,alfa1,alfa1,alfa3)= byxxx(:,alfa1,alfa2,alfa3);
                    byxxx(:,alfa1,alfa3,alfa1) = byxxx(:,alfa1,alfa2,alfa3);
                    byxxx(:,alfa3,alfa1,alfa1) = byxxx(:,alfa1,alfa2,alfa3);                
                end
                % Using symmetry for alfa2 and alfa3            
                if alfa1 ~= alfa2 && alfa2 == alfa3  %alfa1~=alfa2==alfa3 
                    %byxxx(:,alfa1,alfa2,alfa2)= byxxx(:,alfa1,alfa2,alfa3);                
                    byxxx(:,alfa2,alfa1,alfa2) = byxxx(:,alfa1,alfa2,alfa3);                
                    byxxx(:,alfa2,alfa2,alfa1) = byxxx(:,alfa1,alfa2,alfa3);                                
                end            
                % Using symmetry for alfa1,alfa2, and alfa3            
                if alfa1 ~= alfa2 && alfa1 ~= alfa3 &&  alfa2 ~= alfa3 %alfa1~=alfa2~=alfa3
                    %byxxx(:,alfa1,alfa2,alfa3) = byxxx(:,alfa1,alfa2,alfa3);
                    byxxx(:,alfa1,alfa3,alfa2) = byxxx(:,alfa1,alfa2,alfa3);
                    byxxx(:,alfa3,alfa1,alfa2) = byxxx(:,alfa1,alfa2,alfa3);
                    byxxx(:,alfa3,alfa2,alfa1) = byxxx(:,alfa1,alfa2,alfa3);
                    byxxx(:,alfa2,alfa3,alfa1) = byxxx(:,alfa1,alfa2,alfa3); 
                    byxxx(:,alfa2,alfa1,alfa3) = byxxx(:,alfa1,alfa2,alfa3);                 
                end
            end
        end
    end
end

modelYield                    = model;
modelYield.bySS               = bySS;
modelYield.by0                = by0;
modelYield.byx                = byx;
modelYield.byxx               = byxx;
modelYield.byxxx              = byxxx;
modelYield.Byxx               = reshape(modelYield.byxx,maxMat,nx^2);
modelYield.Byxxx              = reshape(modelYield.byxxx,maxMat,nx^3);
modelYield.MAEbonds           = MAEbonds;
modelYield.RMSEbonds          = RMSEbonds;
modelYield.MaxAEbonds         = MaxAEbonds;
modelYield.thetaBondPrice     = thetaBondPrice;
end

